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Both topological crystalline insulators surfaces and graphene host multi-valley massless Dirac 
fermions which are not pinned to a high-symmetry point of the Brillouin zone. Strain couples 
to the low-energy electrons as a time-reversal invariant gauge held, leading to the formation of 
pseudo-Landau levels (PLL). Here we study periodic pseudo-magnetic helds originating from strain 
superlattices. We study the low-energy Dirac PLL spectrum induced by the strain superlattice and 
analyze the effect of various polarized states. Through self-consistent Hartree-Fock calculations we 
establish that, due to the strain superlattice and PLL electronic structure, a valley-ordered state 
spontaneously breaking time-reversal and realizing a quantum Hall phase is favored, while others 
are suppressed. Our analysis applies to both topological crystalline insulators and graphene. 


I. INTRODUCTION 

The discovery of graphene and topological insula¬ 
tors has significantly boosted the ubiquity of condensed 
matter realizations of Dirac fermions as emergent elec¬ 
tronic excitations at low-energy m- Dirac electrons 
in condensed matter systems have enjoyed an enor¬ 
mous amount of interest both from a fundamental and 
technological application perspective [4]. A key differ¬ 
ence between graphene and topological insulators is the 
number of species, or valleys, of Dirac fermions and 
their locations in momentum space. Topological insu¬ 
lators (TI) protected by time reversal symmetry host a 
single-valley Dirac fermion, which is pinned to a time- 
reversal-invariant (TRI) momentum in the surface Bril¬ 
louin zone [5]. In contrast, graphene hosts two valleys of 
Dirac fermions located at non-TRI momenta EHZ], each 
valley having an additional spin degeneracy. More re¬ 
cently, a new type of Dirac fermions was discovered on the 
surface of topological crystalline insulators (TCI) SnTe, 
(Sn,Pb)Se and (Sn,Pb)Te [5HT^. which are protected by 
mirror symmetry of the crystal [HlIIMi]- These Dirac 
fermions exhibit spin-momentum locking as in TIs; how¬ 
ever, there is an even number of Dirac cones at non-TRI 
momenta, a feature similar to graphene. 

In general, when Dirac points are located at non-TRI 
momenta, nonmagnetic perturbations such as strain are 
able to move Dirac points in momentum space, thereby 
acting as an effective gauge held on Dirac fermions. For 
example, strain induces opposite gauge helds for the two 
Dirac valleys in graphene, and spatially inhomogeneous 
strain gives rise to effective magnetic helds that are op¬ 
posite in two valleys, preserving time reversal symme¬ 
try [TTHini- the presence of such a pseudo-magnetic 
held B, the low-energy electronic structure takes the form 
of pseudo-Landau levels (PLLs) with energies charac- 
teric of Dirac electrons in magnetic helds, i.e. ^ \/riB, 
where n is the Landau level (LL) index. Key signatures 
of pseudo-magnetic helds have been experimentally ob¬ 
served in graphene HOI HI]. 

The PLLs have a large single-particle degeneracy, 


which makes them susceptible to many-body instabilities 
in a manner similar to magnetic-held induced LLs [22j . 
Electronic interactions are expected to lift the degener¬ 
acy and drive the system into various gapped states. Two 
primary examples are spin-polarized and valley-polarized 
states of PLLs in graphene [23U28] . 

In this work we consider interacting Dirac electrons 
under periodically modulated pseudo-magnetic helds, 
where regions of positive and negative helds alternate 
in space, forming a superlattice. This held prohle leads 
to a novel electronic structure markedly diherent from 
uniform pseudo-magnetic helds [29l [30]. There are 
various ways in which periodic pseudo-magnetic helds 
can arise, one prominent way being a strain superlat¬ 
tice in graphene or TCI. Such spatially periodic strain 
helds are particularly relevant, as they were experimen¬ 
tally found to develop at interfaces of heterostructures 
built from TCIs (e.g., SnTe) and trivial insulators (e.g., 
PbTe) [3TI - l33| . At these interfaces, the lattice con¬ 
stant mismatch causes dislocations which self-organize 
into a periodic array and therefore produce a natural 
realization of periodic strain helds. The key charac¬ 
teristic of the corresponding periodic pseudo-magnetic 
helds is that they can exist over macroscopic regions. 
In contrast, uniform pseudo-magnetic helds cannot ex¬ 
ist in the thermodynamic limit, owing to the bounded¬ 
ness of strain (=pseudo-gauge held) |30|, unlike a real 
magnetic held. Apart from strain superlattices, periodic 
pseudo-magnetic helds can arise as a result of incom¬ 
mensurate electrostatic potentials originating from, for 
instance, from lattice-mismatched substrates with a twist 
angle [5111551BT] . 

Starting from the strain-induced pseudo-magnetic su¬ 
perlattice, we address the effect of electron-electron in¬ 
teractions. Spatially alternating pseudo-magnetic helds 
change the low-energy electronic structure close to the 
Dirac points. Most strikingly, the energy-momentum dis¬ 
persion in the vicinity of each Dirac point becomes nearly 
hat, leading to a segment of hat band with a twofold de¬ 
generacy. These hat bands arise from the zeroth PLL 
in regions with strong pseudo-magnetic helds; and the 
twofold degeneracy corresponds to Landau orbitals that 
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reside in different spatial regions of opposite fields and 
have opposite Dirac spinor components [30] . a novel fea¬ 
ture that is absent in the case of uniform magnetic fields. 
Counting the two valleys, the flat bands have fourfold 
degeneracy. The presence of flat bands leads to a diverg¬ 
ing density of states, in contrast to the vanishing density 
of states at the Dirac point of massless Dirac fermions. 
Consequently, periodic strain fields provide a feasible and 
effective way of engineering density of states, i.e., elec¬ 
tronic compressibility, at zero energy. 

At charge neutrality, the degenerate flat bands are 
mainly responsible for driving the spontaneous forma¬ 
tion of ordered states. We discuss the various possibilities 
for degeneracy lifting in the flat band and discriminate 
between energetically favorable and unfavorable states. 
Two prominent candidate ordered states are the charge- 
ordered state, where charge is redistributed from the re¬ 
gion of positive (negative) to negative (positive) pseudo- 
magnetic field, and the valley-ordered state, where in 
each spatial region the valley degeneracy is lifted. We 
show that the valley-ordered state in graphene and TCIs 
spontaneously breaks time-reversal symmetry and real¬ 
izes an integer quantum Hall effect similar to the Hal¬ 
dane state [7], but with the important difference of be¬ 
ing driven by electron interactions without any external 
time-reversal-breaking field. 

We determine the mean-field ground state by self- 
consistently solving the full gap equation of interacting 
Dirac fermions under a periodically alternating pseudo- 
magnetic field. The continuum Hamiltonian is micro¬ 
scopically implemented using a lattice model with a 
strain superlattice. Our analysis shows that the support 
of the flat band wavefunctions is of great importance. 
Flat bands in any spatial region only have a twofold val¬ 
ley degeneracy, protected by the time reversal symmetry. 
Therefore lifting this degeneracy by interactions implies 
time reversal symmetry breaking. For this reason, we find 
the valley-ordered quantum Hall state is greatly favored 
over the charged-ordered state under generic forms of 
electron interactions. We show that the order parameters 
corresponding to the ordered states follow the strain pro¬ 
file, highlighting the crucial role of the pseudo-magnetic 
field. 

The present setup for a spontaneous time-reversal sym¬ 
metry breaking quantum Hall state relying on a strain- 
induced flat band should be contrasted with the proposed 
Haldane mass generation for interacting massless Dirac 
fermions in graphene |43H45) . Exact diagonalization and 
DMRG studies [I5H5T| seem to have failed to find the 
interaction-driven Haldane phase in models so far pro¬ 
posed, in contradiction to Hartree-Fock results. It is be¬ 
lieved that the absence of the Haldane phase in ED and 
DMRG phase diagrams stems from the vanishing density 
of states at the Dirac point, and the resulting absence of a 
weak-coupling instability. In contrast, the quantum Hall 
state in our setup already occurs spontaneously at weak 
coupling, owing to the strain-induced flat band. 


II. DIRAC FERMIONS IN TWO DIMENSIONS 

We set out to study the coupling of time-reversal in¬ 
variant pseudo-gauge fields to Dirac electrons. With two 
specific realizations in mind, graphene and TCI surface 
states, we focus on Dirac electrons in two dimensions. 
Let us start by stating the essential features of Dirac 
electrons coupled to pseudo-gauge fields, independent of 
specific context. 

Any two dimensional system respecting time-reversal 
invariance and having Dirac fermions not pinned to a par¬ 
ticular time-reversal invariant momentum, will consist of 
two species of Dirac fermions. Labeling the two species 
by -|- and —, the Dirac Hamiltonian describing the two 
species takes the general form 

ii± = ±Lhv dx — iT^dy)^± ( 1 ) 

where r* is a set of Pauli matrices acting on the pseu¬ 
dospin degree of freedom of the Dirac fermions. Time- 
reversal symmetry relates the two species by exchanging 
4'+ T_. 

Pseudo-gauge fields couple to the Dirac fermions in a 
manner similar to real electromagnetic gauge fields, with 
one crucial difference, however. In order to respect time- 
reversal invariance, the pseudo-gauge field must couple 
to the fermions in such a way that the two species see 
opposite fields. As a result, in the presence of a pseudo¬ 
gauge field given by Ay, {p = x,y), the Dirac Hamiltonian 
becomes 

'H± = ±VF'^±T^{Py ± Ay)'^±. ( 2 ) 

This Hamiltonian (with py = —idy) describes the generic 
Dirac electrons coupled to pseudo-gauge fields. Pseudo- 
Landau level quantization will occur when the gauge field 
Ay acquires spatial dependence, i.e.. Ay = Ay(r). 

The interpretation of the Dirac fermion pseudospin and 
valley degrees of freedom will depend on the particular 
realization of pseudo-magnetic field coupling in a given 
material. In this work we will discuss two examples of 
low-energy Dirac electrons coupled to time-reversal in¬ 
variant gauge fields, which we introduce in the remainder 
of this section. First we consider the case of graphene, 
and then we consider surface states of TCIs. Whereas 
in graphene the Dirac pseudospin degree of freedom de¬ 
rives from the two sublattices [6], the pseudospin of the 
TCI surface state is more complicated due to intrinsic 
spin-orbit coupling, as we will discuss below.Importantly, 
in case of the latter, spin-orbit coupling leads to spin- 
momentum locking in the surface state Dirac theory. 

In both cases, graphene and TCI, the emphasis will be 
on strain-induced pseudo-magnetic field coupling. How¬ 
ever, we will use the case of graphene to point out that 
pseudo-magnetic fields can have a physical origin differ¬ 
ent from strain, giving way to an even wider application 
of our results. 
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A. Dirac fermions in graphene 


The low-energy theory of graphene at charge neutral¬ 
ity is one of the hallmark examples of a 2D Dirac the¬ 
ory [1 m 137]. The two species of nodal Dirac fermions 
are located at the two inequivalent BZ corners, i.e., the 
Dirac points or valleys, and are labeled by K+ and 
K- corresponding to the momenta K+ = (47r/3,0) and 
K_ = — (47r/3,0). The Dirac Hamiltonian is obtained by 
expanding the band structure around the Dirac points 
K and K' in small momenta q relative to the Dirac 
points [6|. It is given by 

'H(g) = + qyT^) = HvFqf^Tf, (3) 


(where vf = V3ta/{2H)). The set of Pauli matrices r* 
acts on the sublattice degree of freedom {A/B) and the 
set of matrices v'^ acts on the valley degree of freedom 
{K+fK- ). In addition, we have defined the Dirac ma¬ 
trices Ta; = and Tj^ = The Hamiltonian acts 

on the Dirac spinor defined by 




■fpA+iq) ■fpB+{q) 'ipA-iq) 



(4) 


Note that we choose the basis so that the A and B lattice 
are exchanged in the K- valley, meaning that we are 
working in the chiral representation (i.e., 'H± = ±(f-T for 
valley K±). 

Starting from the low-energy Dirac Hamiltonian of 
Eq. (§, we introduce a generalized time-reversal invari¬ 
ant pseudo-gauge field by coupling the Dirac fermions to 
the field 


A^ = {A^,,Al), (5) 

which consists of three components i = 1,2,3. The cou¬ 
pling to the fermions has the same form as ordinary min¬ 
imal coupling, but with different gauge charges D® ex¬ 
pressed as 

n{q) = hvFr^iq^ + A^^n^). ( 6 ) 

The gauge charge matrices fl® encode the distinct nature 
of the pseudo-gauge field as compared to the ordinary 
gauge field, and are given by H® = The 

third gauge charge matrix = A is diagonal in valley 
space and assigns opposite sign to the two valleys. There¬ 
fore, the component A^^pP realizes the general Hamilto¬ 
nian of Eq. in graphene. In graphene, this is the 
pseudo-gauge field component that arises in the presence 
of strain and plays a central role in this work. The pres¬ 
ence of a field A coupling to leads to a moving of the 
Dirac points away from and Ar_, in opposite direc¬ 
tions. This is shown in Fig. (left), where the bold blue 
dots denote the Dirac points moving towards the zone 
center. 

The following properties of the gauge field charges will 
be important for our analysis. The charges H® real¬ 
ize a pseudospin SU{2) algebra, expressed as = 


2iP^Vl.^. The matrices D® commute with the Hamilto¬ 
nian in the absence of fields, and as a consequence gen¬ 
erate a continuous SU(2) symmetry of the low-energy 
graphene Hamiltonian. This symmetry is broken when 
mass terms are introduced to the Hamiltonian, i.e., when 
the Dirac electrons are gapped out. In particular, the 
set of mass matrices T = (ri,r 2 ,r 3 ) = 
describes masses that anti-commute with the Hamilto¬ 
nian and between themselves. They constitute a set of 
compatible masses, the physical nature of which is well- 
known. Specifically, the mass Ts = corresponds 

to an electrostatic potential making the two honeycomb 
sublattices inequivalent and breaking inversion symme¬ 
try. Such term is diagonal in valley-space, i.e., it does not 
couple the two Dirac points. The other two masses, Ti 
and r 2 , which are off-diagonal in valley space, are known 
as Kekule masses and correspond to modulations of the 
tight-binding nearest-neighbor hopping parameter t with 
tripled unit cell [31]. The breaking of translational in¬ 
variance and the modulations over small distances (large 
momenta) couple the Dirac points. 

The gauge charges H® act as generators of rotations 
within the space of masses, which follows from the com¬ 
mutation relation [fl®,rj] = 2ie^^’^Tk. In addition to the 
mass terms F, there is a mass term the time-reversal 
odd Haldane mass |7], which anti-commutes with the 
Hamiltonian <§ , but commutes with both the F^ and 
the H®. Hence, whereas F is a vector under the transfor¬ 
mations generated by H®, is a scalar. 

The two remaining gauge charges and 

are off-diagonal in valley space but diagonal in 
sublattice space. The former implies translational sym¬ 
metry breaking and the latter implies that these terms 
arise due to charge density modulations. Consequently, 
charge density waves (CDWs) with a six-site unit cell, 
which we will refer to as valley-coupling CDWs, lead to 
a pseudo-gauge coupling in the same way as strain [34]. 
The SU{2) structure of the gauge charges H® implies that 
within the low-energy theory, the pseudo-gauge field com¬ 
ponents are unitarily equivalent to each other. 


B. TCI surface state Dirac fermions 

Topological insulator materials are bulk insulators 
hosting gapless Dirac fermions at their surfaces mm- 
The spin-momentum locked surface Dirac fermions are 
protected by time-reversal symmetry, and as a result 
they are pinned at the time-reversal invariant momenta 
(TRIM). Due to this symmetry-protected pinning, the 
surface states of topological insulators do not allow for 
time-reversal invariant pseudo-gauge field coupling. In 
particular, strain is not able to move the Kramer’s dou¬ 
blet away from the TRIM. 

In contrast, the TCIs are topological materials pro¬ 
tected by crystalline symmetries mm, which host sur¬ 
face Dirac fermions not pinned to the TRIM [TS] |3S1 [SS] ■ 
As a result, strain can couple to the low-energy Dirac 
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FIG. 1. (Left) Hexagonal Brillouin zone of the honeycomb 
lattice. The two Dirac points K and K are marked by bold 
blue dots. The blue arrows indicate possible Dirac points 
moving towards the zone center due to strain, i.e. ~ Uxx — Uyy. 
(Right) Square surface Brillouin of TCI surface state with 
two sets of Dirac points located at Xi (blue) and X 2 (red). 
Arrows indicate the moving of Dirac points towards the zone 
center due to symmetric strain ~ Uxx + Uyy. 


fermions as a pseudo-gauge field and can move the Dirac 
points in momentum space, in a way that depends on the 
symmetry of the strain tensor [niisn]- 

In this work we specfically focus on the SnTe mate¬ 
rial class [8] and its mirror symmetry-protected surface 
Dirac fermions appearing on the (001) surface. The sur¬ 
face Brillouin zone of the (001) surface is shown in Fig.[^ 
Two species of low-energy Dirac fermions related by time- 
reversal symmetry exist in the vicinity of the surface 
time-reversal invariant momenta Xi and ^ 2 , represented 
as blue and red dots in Fig. The surface state Dirac 
Hamiltonian at Xi, given by the terms that respect the 
crystal symmetries leaving Xi invariant, reads 

TLx^ {q} = viqia"^ — V 2 q 20 '^ + mv^ + (7) 

and a similar expression can be derived for AI 2 . Here a’’ is 
a set of Pauli matrices that represents a Kramers doublet, 
and zz* is a valley degree of freedom corresponding to the 
two inequivalent bulk L-points mapped onto Xi. The 
momentum q is measured with respect to Xi] the spin- 
momentum locking shown in Hamiltonian Q (i.e., first 
two terms) comes from spin-orbit coupling. For m = 
(5 = 0 there are two degenerate Kramer’s doublets at 
Xi, which are split in energy by finite m and 6. Most 
importantly, finite m and 5 leads to the appearance of 
two species low-energy Dirac points, which are located 
at A± = (0, + (5^/u 2), measured from Xi. 

The Hamiltonian of Eq. can be projected into the 
subspace corresponding to A± to obtain the effective low- 
energy Dirac theory. This yields m 

(9) = ( 8 ) 


tonian is valley-isotropic, taking iz’” = ±1 as an effective 
valley degree of freedom. 

With the Hamiltonian of Eq. ([^ we have arrived at 
a desciption of the low-energy Dirac fermions that has 
the general form introduced in the beginning of this sec¬ 
tion, and is thus similar to the graphene Hamiltonian of 
Eq. (§. Hence, in a way analogous to graphene, we can 
use symmetry arguments to establish the effect of various 
perturbations. For instance, since the TCI surface states 
are protected by mirror symmetry, one expects mirror 
symmetry breaking to open up a gap. We find two such 
gap opening mass terms, which do not couple the low- 
energy valleys, and they are given in the basis of (|^ as 

and jz^T^. The former is a time-reversal even mass and 
corresponds to a ferroelectric distortion of the crystal. It 
derives from the Dirac bilinear in the basis of ([^. The 
mass term zz^r^ breaks time-reversal symmetry and orig¬ 
inates from the terms and cr^zz^ in the basis of Eq. Q. 
The mass gap originating from zz^r^ is equivalent to the 
graphene Haldane gap, and consequently corresponds to 
a QAH phase [IS]. 

Similarly, by using symmetry arguments, the time- 
reversal invariant pseudo-gauge field couplings can be 
identified. As a consequence of the low symmetry of 
the Xi point, there are no two dimensional representa¬ 
tions which directly imply pseudo-gauge coupling. How¬ 
ever, since the symmetric terms and cr’^zz^ displace 
the Dirac points in momentum space, any perturbation 
coupling to them, will have the effect of a pseudo-gauge 
field. Looking for other terms both even under time- 
reversal and inversion (as expected for strain), one finds 
another Dirac bilinear given by cr^zz^. We will show in the 
next section, when we discuss strain and strain superlat¬ 
tices, that components of the strain tensor couple to these 
terms. The effect of these terms is shown schematically 
in Fig.(right), where bold blue and red dots denote the 
Dirac points in the vicinity of Xi and X 2 , respectively, 
shifting towards the zone center as a result of strain. 


III. PERIODIC STRAIN SUPERLATTICES 


In the previous section we introduced pseudo-gauge 
field coupling in the context of graphene and TCI sur¬ 
face states, and argued that strain realizes such coupling. 
We now turn to a more detailed discussion of strain, and 
more specifically periodic pseudo-magnetic fields arising 
due to periodically varying strain, i.e., a strain superlat¬ 
tice. 

Elastic deformations of the crystal lattice, i.e., strain 
fields, are described by the strain tensor Uij given by 

Uij — “t“ (9) 


where t* is the effective pseudospin degree of freedom, 
q is now measured with respect to A±, and v{ = 
Vi5j\/m? + (5^. Note that in the chosen basis the Hamil¬ 


where Ui {i = x, y) is the displacement field. Given the 
symmetry of the crystal lattice, the strain tensor can 
be decomposed into components transforming as distinct 
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representations of the symmetry group. From this de¬ 
composition one can read off which lattice deformations 
couple as (pseudo-)gauge fields to the Dirac fermions. 

In case of hexagonal symmetry, applicable to graphene, 
there are two d-wave strain components {uxx — 
Uyy,—2uxy) ~ (dx^-y^^dxy), which couple to the Dirac 
fermions as the valley-diagonal field of Eq. § (let us 
omit the label 3, i.e. A = A^ ^ for the moment) |40j. Note 
that assuming full hexagonal symmetry implies a single 
coupling constant, 

C Uxx — Uyy \ 

^lUxy J 

For illustration purposes, the effect of finite and constant 
Ax is shown graphically in Fig. (left), where the Dirac 
points iF_|_ and K_ move along the kx axis. In case of 
the square symmetry, which applies to the (001) surface 
states of TCIs, the d-wave components dx 2 _x'^ an dxiX 2 
are not degenerate, and one finds at Xi m, 

f f Oii(^Uxx ^yy) 

J Y 0^2Uxy 

In addition, in the previous section we observed that 
a perturbation respecting all symmetries can move the 
Dirac points in momentum space, implying that Uxx+Uyy 
enters the expression for Ax as well, with an independent 
coupling. It is the latter strain component, Uxx + Uyy, the 
effect of which is shown in 0 (right). 

Let us now come to the case of periodic strain fields 
with wavelength A. More specifically, we consider Uij —t 
Uij{r) and consequently A — >■ A{r). The periodicity of 
A{r) is directly reflected in the periodicity of the pseudo- 
magnetic field B = V X A(r) = B{r), which should 
be compared to and contrasted with uniform pseudo- 
magnetic field B. In order to implement the strain super¬ 
lattice in a tight-binding setting, we take the graphene 
lattice as a simple regularization of the continuum the¬ 
ory. To solve the superlattice Hamiltonian, we estab¬ 
lish a connection between the strain components and the 
change in overlap intergals Stn, where n = 1,2,3 labels 
the three nearest neighbor vectors {i5„}. The overlap in¬ 
tegral change is expressed in terms of the strain tensor 
Uij as 6tn = dhdi,Uij, which becomes 

'^xx ^yy\ j 2dti 5t2 ^^3 

-2uxy J y \/3{St2 - Sts) 

This expresses the pseudo-gauge field in terms of the 
modulation of the hopping —>■ t -|- 6tn. 

We proceed to consider a single-propagation vector 
pseudo-gauge field modulation and obtain the electronic 
spectrum. A particularly convenient choice is the prop¬ 
agation vector 6G = G/A, where G = (0,47r/-\/3) is a 
reciprocal lattice vector. Then A is the superlattice wave 





length, given in terms of graphene unit cells, i.e. A = 700 
leads to a superlattice unit cell containing 700 graphene 
unit cells. The pseudo-gauge field A and corresponding 
pseudo-magnetic field are given by 

Ax{f) = Acos{6M ■ f) 

B{-F) = -dyAx = SMyAsm{6M ■ r), (13) 

while Ay = 0, since we are only interested in the trans¬ 
verse component. A denotes the amplitude of the strain 
field, i.e., the maximal change in overlap 5t. 

The spectra in the presence of the strain superlattice 
for some values of A are shown in Fig.[^ We observe that 
for increasing A, implying increasing pseudo-magnetic 
field strength, a flat zero-energy band forms at the Dirac 
points, in addition to higher energy dispersive but doubly 
degenerate bands. This specific reorganization of the low- 
energy electronic spectrum resembles the Landau level 
structure of external magnetic fields. We will establish 
a detailed connection between Landau level physics and 
periodic strain in the next section. A key feature we wish 
to stress here, is that the formation of the zero-energy flat 
band, the degeneracy of which is related to the strength 
of the pseudo-magnetic field, leads to a finite and consid¬ 
erable density of states (DOS) at the charge neutrality 
point. In stark constrast, in the unstrained case Dirac 
electrons have linearly vanishing DOS at the charge neu¬ 
trality point. 

Instead of the propagation vector JG, we can take the 
propagation vector 5K = AT/A, where K is the Dirac 
point vector defined above. This may be viewed as 
a simple rotation of 5G, which will result in modified 
Moire pattern. We then have for the spatially dependent 
pseudo-gauge field 

Ax{'G) = Acos{SK ■ r) 

B{r) = —dyAx = 5KyAsm{5K ■ r) (14) 

where is it important to choose K such that the field 
has a nonzero transverse component. For given A the 
strain superlattice unit cell constains 3A graphene unit 
cells, which has the benefit that it is commensurate with 
any perturbation modulated by K coupling the Dirac 
points. This allows to treat valley-diagonal and valley- 
off diagonal perturbations on the same footing. 

We recall that the low-energy Dirac Hamiltonian in the 
presence of strain reads % = hvpB y{—idy+A?y^{r)BL^). As 
noted above, a unitary matrix U can be used to rotate to 
another gauge field component, WHU = hvpB+ 
Al^(r)f2^). Clearly, this does not change the spectrum 
and therefore electrostatic potential superlattices, which 
would couple to the gauge field components and 
[n] , are equivalent to strain superlattices. Therefore, 
even though we focus on strain in this work, we high¬ 
light that in the context of graphene spatially modulated 
valley-coupling CDWs induce periodic time-reversal in¬ 
variant pseudo-magnetic fields in the same way as strain. 
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FIG. 2. Spectra of graphene in the presence of periodically modulated strain for different values of the amplitude of modulation 
(given by A ~ 25ti — 5t2 — see main text). The modulation wavelength A is chosen to be 700 graphene unit cells and the 
propagation vector is 5G = (0, 47r/\/3)/A. The amplitudes are (a) A = O.Olf, (b) A = 0.03t, (c) A = 0.05t. Note that the plots 
show K' folded onto the fc^j-axis. 


IV. FLAT BAND PSEUDO-LANDAU LEVELS 
INDUCED BY STRAIN 


and the dynamical momenta obey the commutation re¬ 
lations 


In this section we address the spectral properties of 
Dirac electrons in the presence of a time-reversal invari¬ 
ant pseudo-magnetic field. As in the previous section, 
we particularize to the case of graphene (i.e., our lat¬ 
tice regularization of the continuum theory), and start 
by considering a spatially uniform field induced by strain. 
Uniform fields are fundamentally different from periodi¬ 
cally modulated fields induced by the strain superlattice, 
but we can use the results for the former to develop an 
intuition for the case of periodic pseudo-magnetic fields. 
In particular, we may, for the sake of argument, think of 
the periodic field as alternating regions of positive and 
negative constant fields, which is schematically shown in 
Fig. 15 (top). 


A. Uniform pseudo-magnetic fields and PLLs 


The Hamiltonian is given by Eq. j for which we only 
retain the strain component A = A^, 

= hvF^xiqx + Ax^^) + hvpTyiqy + AyV,^), (15) 


and A = A{x) is taken so as to describe a constant field. 
The mathematical structure of the Hamiltonian is equiv¬ 
alent to that of a real eletromagnetic field and we can 
use known techniques to solve it. For completeness we 
briefly review the essentials here, leaving more details to 
Appendix]^ We first introduce dynamical momenta H^ 
and H^ for each valley iz = ±, reflecting the fact the sign 
of the pseudo-magnetic field is opposite for the valleys 
(recall that = v^). The Hamiltonian for each of the 
valleys reads 



[n±,n±] = T*Jf. (17) 

These commutation relations can be used to define rais¬ 
ing and lowering operators in the standard way (see 
App.|g , in terms of which the Hamiltonian takes the 
form 


n± = ± 




(18) 


Here we have defined and the the rais¬ 

ing and lowering operators obey the commutation rela¬ 
tion [a±,a^] = ±1. This commutation relation is a key 
feature of time-reversal invariant pseudo-magnetic fields, 
since it reflects anti-parallel field alignment in the two 
valleys. The operation of raising and lowering is inter¬ 
changed for the two valleys, which has important conse¬ 
quences for the structure of the eigenstates. In particular, 
the eigenstates of the PLL zero modes are localized on 
the same sublattice, instead of on opposite sublattices. 
More specifically one finds 

We stress that this implies localization on the same sub¬ 
lattice, given the choice of basis in Eq. (j^). Time-reversal 
symmetry is preserved by counterpropagation of Landau 
orbitals in the two valleys (i.e., |</5g j.) = |</5 q _j.)). The 
n = 0 PLL has energy E = 0. Eigenstates corresponding 
to n ^ 0 PLLs take the form 


14'+ 


n±/ 


1 / \^n—l,k) 

i ±|v5n,fe) 




^ ( \<,k) \ 

\T\V’n-l,k)) 

( 20 ) 


'H±{q) = ±vf 


- itiy 


(16) 
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FIG. 3. Upper left: Spectrum in the presence of a 
strain-induced periodic pseudo-magnetic field, shown with 
both Dirac points folded onto F. Periodicity of the pseudo- 
magnetic field is 1200 graphene unit cells, and we used A = 
O.OSt. Black arrows indicate for which k the zero mode eigen¬ 
states are plotted in the upper and lower right panels. In these 
two panels we show the wave function distribution ~ k\'^ 

of the full zero mode (i.e. n = 0) subspace over the graphene 
unit cells for the A-sublattice (lower right, red) and B- 
sublattice (upper right, black). Black arrows indicate which 
k they correspond to in the upper left plot; the blue arrow in¬ 
dicates in which order. Lower left: plot of the periodic strain 
modulation Ax ~ Acos(27rj//A) ~ 2Sti{y) — 5t2{y) — 5t3{y) 
(red) and corresponding pseudo-magnetic held. Note that for 
clarity we have rescaled the amplitude of the pseudo-magnetic 
held to A. 
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and they have energies E±{n) = ± 1 / 2 ^^ for each valley. 

In order to gain insight into the effect of perturba¬ 
tions on periodic strain superlattices, we review the ef¬ 
fect of such perturbations on the PLL spectra for uniform 
pseudo-magnetic field we just presented. Consider the 
n = 0 PLL, the eigenstates of which are given in Eq. (19). 
Let us first comment on the mass gap terms and 

T*. The effect of the sublattice CDW and the TRS 
breaking Haldane term is reversed as compared to real 
magnetic fields [23]. The sublattce polarized CDW sim¬ 
ply shifts the energy of zero modes but does not break 
their degeneracy, whereas the Haldane mass energet¬ 
ically splits the zero modes in a symmetric way. This 
is an immediate consequence of the sublattice structure 
of the PLL zero modes. The two Kekule masses and 
do not affect the zero modes at all, they are neither 
split nor shifted, since they are off-diagonal is sublattice 
space. 

Perturbations that do split the zero modes in a fash¬ 
ion similar to the Haldane term are charge density waves 
with tripled unit cell, i.e., charge density waves that cou¬ 
ple the valleys iL+ and [Mj. These charge density 
waves couple to the Dirac mastrices ^ iP't^ and 

iP'T^. Projecting these into the PLL zero mode subspace, 
one finds effective Pauli matrices and f^, which anti¬ 
commute with the Haldane mass projected into the zero 


mode space, —)• f^. This leads to the counterintuitive 
situation of anticommuting masses only one of which is 
TRS breaking and nontrivial |32|. We note that these 
charge density waves with tripled unit cell correspond 
to the other gauge field components of the SU{2) gauge 
field, i.e. they enter as A^., A^, A^, and Ay in Eq. §. 
Yet another perturbation that splits the zeroth PLL is 
the valley mass, given by A ^ making the valleys inequiv¬ 
alent, but acting as the identity in sublattice space. Its 
spectral effect is equivalent to that of the Haldane term, 
meaning a symmetric splitting of the zero modes. 

Understanding the spectral effect of these Dirac 
fermion bilinears on the zeroth PLL gives a first idea 
of the ways in which their spontaneous formation can 
lower the energy for charge neutral systems. For a more 
refined understanding of the energetics it is necessary to 
consider the effect of perturbations on higher PLLs (i.e., 
n 7 ^ 0). For both the CDW term and the Hal¬ 
dane term all PLLs with n 7 ^ 0 get pushed up or 
down in energy depending on whether they have energies 
±-\/ 2 ^^n, i.e. positive (negative) solutions get pushed up 
(down). This is different for the valley mass A, which 
pushes all PLLs of valley K+ up and of valley K- down, 
effectively splitting all PLLs, even the n 7 ^ 0 levels, leav¬ 
ing do degeneracies behind. The charge density waves 
with enlarged unit cell (i.e., coupling the valleys) both 
split and shift the higher order PLLs, which may be seen 
straightforwardly by using perturbation theory up to sec¬ 
ond order. Based on these considerations we obtain an 
intuition for the spontaneous generation of Dirac fermion 
bilinears due to interactions, depending on the location 
of the Fermi level. 


B. Alternating pseudo-magnetic fields and 
superlattice PLLs 


The assumption of uniform pseudo-magnetic held is a 
useful hrst step towards understanding the physics of su¬ 
perlattice PLLs. As a next step we consider the case 
of alternating pseudo-magnetic held, which for simplic¬ 
ity we will take as a periodic arrangement of regions of 
positive and negative constant held. This will provide 
valuable insight in the case of periodic harmonic pseudo- 
magnetic helds. Schematically this held arrangement is 
shown in the top row of Fig. Figure [^ shows how we 
can think of the a periodically alternating held as a strain 
superlattice with ehective “two-site” unit cell (i.e., posi¬ 
tive and negative held), reminiscent of an antiferromag- 
net, leading to a doubling of the PLL degeneracies. For 
instance, the space of zero mode PLLs is doubled, since 
we have the spatial degeneracy in addition to valley de¬ 
generacy. For each valley there is a zero mode localized 
in the positive held region, meaning on the A sublattice, 
and a zero mode localized in the negative held region, on 
the B sublattice. 

The additional degree of freedom originating from the 
periodicity of the pseudo-magnetic held gives rise to a 
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FIG. 4. Schematic representation of the PLL structure in 
case of spatially alternating pseudo-magnetic field. The al¬ 
ternating regions of positive and negative field are depicted 
in the upper row. In the row below we depict the degener¬ 
ate n = 0 PLLs in the two regions, the structure of which 
periodically alternates in accordance with the field. Here A 
and B label a general pseudo-spin degree of freedom, which 
corresponds to the A/B sublattice in graphene. The bottom 
two rows show the two prime polarized state candidates. In 
the charge ordered state with ferroelectric polarization, the 
energy levels in the positive field (A) have higher energy than 
the negative held (B) states, the latter being fully occupied. 
In the valley-ordered quantum Hall state (very bottom) lev¬ 
els are split in each region, occupying a single valley in each 
region. 


richer structure of polarized or ordered states. Focus¬ 
ing on the PLL zero mode subspace, relevant at charge 
neutrality, there are multiple ordered states that lift the 
degeneracy of the zero mode subspace. Two of them are 
shown in Fig.|^ The first is a charge ordered state, where 
the zero mode PLLs in one of the two spatial regions are 
both occupied, leaving the zero modes in the other re¬ 
gion unoccupied. This leads to a redistribution of charge 
between the two regions and an associated ferroelectric 
polarization along the propagation direction of the su¬ 
perlattice wave vector. In case of graphene this state 
is realized by the sublattice CDW, which energetically 
discriminates the sublattices. The other state shown in 
Fig. is the valley-ordered quantum Hall (or Haldane) 
state. In such state the zero modes corresponding to 
an “up” pseudo-magnetic field are occupied. Note that 
this implies an alternating occupation of valleys K±^ as 
shown in Fig. Therefore, this state may be called 
anti-ferro-valley-ordered. In graphene the valley-ordered 
quantum Hall state is realized by the time-reversal sym¬ 
metry breaking Haldane term. A third PLL polarized 
state is obtained by occupying the same valley in each 
spatial region. This state also breaks time-reversal, but 
contrary to the valley-ordered quantum Hall state the 
pseudo-magnetic field seen by the occupied PLLs alter¬ 


nates. The inversion of PLL occupation in one of the two 
regions with respect to the anti-ferro-valley-ordered state 
suggest the name ferro-valley-oideied. In graphene such 
state is realized by the valley mass term. 

We now establish a connection between the simplified 
description of alternating pseudo-magnetic fields in terms 
of continuum PLLs, and the periodic strain (super)lattice 
model introduced in the previous section. In order to 
do so we take the unidirectional periodic strain profile 
compatible with an enlarged six-site unit cell (i.e., Dirac 
valleys folded onto F) defined in Eq. (I4|. Solving the 
tight-binding Hamiltonian in the presence of the strain 
superlattice yields the spectrum shown in Fig. (upper 
left). The connection is made by interpreting this spec¬ 
trum in terms of PLLs. 


Let us first study the wave function support of the 
zero energy solutions and compare that to the zeroth 
PLL. The right column of Fig. shows the wave func¬ 
tion support feP of the wave functions ipn=o,k cor¬ 

responding to zero energy solutions (n = 0) labeled by 
k {k should be identified with ky). Black arrows ex¬ 
plicitly indicate which k corresponds to which 
profile. Since there are two valleys and the strain su¬ 
perlattice unit cell consists of two distinct regions with 
opposite pseudo-magnetic field, there is a fourfold de¬ 
generacy at each k. This is reflected in Fig. where two 
pseudo-magnetic Landau-like orbitals are localized on the 
A-sublattice (lower right) and two localized on the B- 
sublattice. The wave function support clearly shows the 
spatial separation of solutions living on distinct sublat¬ 
tices. In the region of positive field (see Fig. [flower left) 
“zero modes” are localized on the A-sublattice, and on 
the H-sublattice in the region of negative field. In addi¬ 
tion, we observe that for k moving away from F (following 
the blue arrow) the two Landau orbitals in each region 
move away from the position of maximum field towards 
the position of vanishing field. In particular, they move 
away in opposite directions, which is a direct consequence 
of their different valley index. As the two valleys effective 
see opposite fields, and the spatial position of a Landau 
orbital is given by a: cx sgn(H)fcy, the Landau orbitals are 
expected to spatially move in opposite directions. With 
increasing k (~ ky) Landau orbitals of different spatial 
regions and the same valley start to overlap, eventually 
leading to the splitting observed in spectrum indicated 
by the most right black arrow in the upper left panel of 
Fig.0 Hence, at the junctions between regions of pos¬ 
itive and negative field, the Landau orbitals acquire a 
dispersion and form a series of snake states [30] • 

Next, we ask what the spectral effect is of the perturba¬ 
tions that are expected to split degeneracies, in particular 
the flat band degeneracies as depicted in Fig. To this 
end, we solve the strain superlattice Hamiltonian in the 
presence of various perturbations, which for the moment 
we take to be spatially uniform, i.e., not follow the super¬ 
lattice envelope, in accordance with the schematic picture 
of alternating regions of constant field (Fig.|^. Figure]^ 
shows the low-energy spectrum in the vicinity of F. The 
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FIG. 5. Spectra of graphene in the presence of periodic strain and in the presence of additional pertnrbations splitting or lifting 
PLL energies. Spectra are obtained for strain snperlattice nnit cells containing A = 600 graphene unit cells and A — O.OSt. (a) 
Free graphene (b) CDW (c) Haldane term (d) valley mass term (e) CDWi with tripled unit cell (f) CDW 2 with tripled unit 
cell. 


unperturbed case corresponding to Hg + .^strain is shown 
in Fig. [^a) for reference. The sublattice charge ordered 
state and (anti-ferro-))valley-ordered quantum Hall state 
are shown in Fig. [^b) and Fig. [^c), respectively. Both 
open up a full gap, splitting the zero mode subspace and 
shifting the higher PLLs, as expected. Figure [^d) shows 
the effect of a valley mass term, i.e., the ferro-valley- 
ordered, on the low-energy spectrum. Whereas in each 
valley degeneracies are preserved, the valleys are split, as 
expected. As a consequence, the spectrum is not gapped 
and the Fermi level crosses the propagating snake states 
associated with the flat band PLLs [3^ . 

Figure [^e) and ^f) show the spectrum obtained in 
the presence of valley-coupling CDWs, which we have 
discussed split the continuum n = 0 PLL in a way sim¬ 
ilar to the Haldane term. We observe that in case of a 
strain superlattice, the CDWs do not lead to a full gap, 
but only lift the degeneracy of the flat band states local¬ 
ized at the position where the pseudo-magnetic field has 
its extrema. The degeneracy is not lifted in the vicin¬ 
ity of the nodes of the periodic pseudo-magnetic held. 
The absence of a full gap in case of periodic strain can 
be understood by considering the spectral effect of the 
valley-coupling CDWs in case of zero pseudo-magnetic 
held. The valley-coupling CDWs do not open up a gap 
in that case, but only shift the Dirac points. Hence, at 
the nodes of the pseudo-magnetic held one expects the 
absence of a gap. Note also that the spectrum numeri¬ 
cally obtained numerically shows both split and shifted 
higher (n 7 ^ 0) PLLs. 


From this analysis we conclude that in the presence of 
the strain superlattice, the low-energy electronic struc¬ 
ture can be approximated by sets of PLLs for the two 
distinct regions of the superlattice unit cell. In addi¬ 
tion, based on their ehect on the degenerate low-energy 
PLLs, we expect the charge-ordered and anti-ferro-valley- 
ordered states to be the dominant instabilities in the pres¬ 
ence of the pseudo-magnetic held superlattice. 

We end this section with two remarks. First, we note 
that the analysis presented here is based on the the as¬ 
sumption of a strain-induced pseudo-magnetic held, i.e., 
A = A^ in Eq. As mentioned in Section]^ a unitary 
matrix can rotate to another component, e.g. A^ or A^. 
This is does not change the (low-energy) spectrum, but 
it does change the nature of the eigenstates. In addition, 
it also changes the nature of the perturbations, since the 
unitary matrix rotates within the space of masses rep¬ 
resented by F as well. In particular, this implies that a 
sublattice polarized term will be rotated into one 
of the Kekule terms. Interestingly, the Haldane term is 
a scalar under these unitary rotations and hence is in¬ 
variant. To summarize, the analysis of this section still 
applies, but in a rotated basis. 

The second remark concerns the applicability of our 
analysis to TCI surface states. The arguments put for¬ 
ward in the present section build on the specific example 
of graphene PLL physics. They remain valid in the con¬ 
text of TCI surface states. Most importantly, in the pres¬ 
ence of a uniform pseudo-magnetic field, the TCI surface 
n = 0 PLLs are localized on the TCI pseudospin de- 
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gree of freedom in such a way that the time-reversal in¬ 
variant ferroelectric distortion of the crystal lattice only 
shifts them, whereas a time-reversal breaking Zeeman- 
type spin-coupling splits them in energy. As a result, 
there is a one-to-one correspondence between the valley- 
ordered quantum Hall state and charge-ordered state in 
graphene, and Zeeman term and ferroelectric distortion 
of TCI surface states. 


V. INTERACTING ELECTRONS IN A STRAIN 
SUPERLATTICE 


In order to systematically investigate the patterns of 
symmetry breaking and PLL splitting, as a consequence 
of interacting flat band electron, we have studied an in¬ 
teracting Hamiltonian on the graphene honeycomb lat¬ 
tice and performed extensive self-consistent Hartree Fock 
calculations. We report the results in this section. 

Based on the intuitive picture of the PLLs in the two 
spatially separated regions we anticipate both the forma¬ 
tion of a charge-ordered state with ferroelectric polariza¬ 
tion, corresponding to a redistribution of charge between 
the regions of positive and negative pseudo-magnetic 
field, and the formation of a valley-ordered quantum Hall 
ground state. Interactions which may lead to the forma¬ 
tion of these states in a graphene lattice model are the 
nearest neighbor (NN) and next-nearest neighbor (NNN) 
density-density interactions [43], respectively, as will be 
demonstrated below. We therefore consider the interact¬ 
ing Hamiltonian H = Hq -|- Hmt where the inter¬ 

acting part Hint is given by 

Hint = ^ flrhr’ + ^2 ^ flrflr'■ (21) 

{rr') 


Here Vi is the NN interaction strength and V 2 the NNN 
interaction strength, and fir is the number operator at 
site r. The sums over (rr') and ((rr')) are over NN and 
NNN, respectively. 

The spatially modulated strain is implemented in the 
way described above, meaning that we take H^ to con¬ 
tain A = (Ax, 0 ), where the Ax component of the strain- 
induced gauge field originates from hopping amplitude 
modulations Ax ~ 25ti — 6 t 2 — which are given a 
spatial dependence dtn —>■ A schematic repre¬ 

sentation of the way we set up the calculations is given 
in Fig. 1^ which is a generalization of the approach of 
Ref. HD to the case of strain superlattices (more details 
in Appendix H- In order to allow for charge density 
waves that couple the Dirac points and K- we work 
with a tripled ( 6 -site) unit cell (red dashed hexagons 
in Fig. [ 6 ) [Hj. The corresponding lattice vectors are 
hi = 23,2 + Qi and 62 = Oi — 02 in terms of the graphene 
lattice vectors Hi. We take the strain-induced gauge field 
to be 


Ax{r) = Acos( 


K' ■ r 
A/3 ’ 



FIG. 6 . Upper left: Schematic representation of the graphene 
Brillouin zone (enter black hexagon) and the folded Bril- 
louin zone corresponding to the tripled unit cell (inner black 
hexagon). Dirac points of pristine graphene are located at 
K and K', which are folded onto F when tripling the unit 
cell. The vectors connecting F to A and K', indicated by 
blue arrows, are reciprocal lattice vectors of the enlarged lat¬ 
tice vectors. Lower left: Wigner-Seitz sells containing three 
elementary graphene unit cells, with lattice vectors ai and 
a 2 . Right: graphene lattice (black hexagons) and tripled unit 
cells (red dashed hexagons), including labeling of sites. Ai 
and Bi label the six sites within the enlarged unit cell; n and 
n-|-l label the position in the strain-induced superlattice. 5tn 
are the hopping amplitude changes due to strain. 


where r denotes the position of an elementary graphene 
unit cell, K is the A_ wave vector as shown in Fig. 
and A equals the number of graphene unit cells in the 
periodic strain superlattice unit cell. The number of 6 - 
site unit cells in the superlattice unit cell is therefore A/3. 
The lattice vectors of the superlattice are hi and A 62 - 
We stress that we consider a fully periodic system with 
wavelength ^ A set by the periodic strain superlattice. 

Spinless electrons on the graphene honeycomb lattice 
serve as a model system for strain-induced PLLs. We 
therefore first present Hartree-Fock results for spinless 
electrons and then briefly comment on the spin degree of 
freedom, recalling that the case of TCI surface states is 
similar to spinless graphene electrons. 


A. Spinless electrons 

Within the framework of standard Hartree-Fock theory 
we decouple the interactions both in diagonal (Hartree) 
and off-diagonal (Fock) channels. The diagonal (or 
charge density) order parameter at a given site in the 
superlattice unit cell is labeled by (a, i, 1), and is defined 
as 


Pai (0 


k 


( 22 ) 


(23) 
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where a = A, B; i = 1,2,3; I labels the 6 -site cell in 
the superlattice cell, and N is the total number of super¬ 
lattice unit cells. We define Pa{l) as the average over i 
within cell I, 


= (24) 

Based on previous considerations we are interested in 
possible local sublattice imbalances Ap_(^) expressed as 

Ap-(0 = Pa(0 - (25) 

and possible ferroelectric redistribution of charges be¬ 
tween regions of positive and negative pseudomagnetic 
field, which can be expressed as 

Ap +(0 = Pa( 0 +Pb (0 - 1 - (26) 

In addition to the charge density order parameter, 
in order to detect the valley-ordered Haldane state, we 
study the quantum Hall (QH) order parameter originat¬ 
ing from spontaneously generated next-nearest neighbor 
tunneling with complex amplitude. Decoupling in the 
off-diagonal Fock channel leads to 

Xhi' = ^ E QXii'ikMS, mcjil', k)), (27) 

k 

where the phase factors Qujii{k) are defined in Ap¬ 
pendix 1^ and from these we extract the QH order pa¬ 
rameter Aqh(0 as 

^Qh (0 = E (28) 

zj7'GNNN(0 


The sum is over all NNN pairs which belong to cell I and 
r]A = —PB = 1 since in the QH state fluxes on the A and 
B sublattices have opposite sign. 

The results of self-consistent HF calculations we 
present here were performed for spinless electrons on lat¬ 
tices of size TV = 16 X 16 and A = 3 x 40 = 120. In the 
following we map out the phase diagram as function of Vi 
and V 2 by discussing the results for three specific regimes 
separately. First we will focus on (Hi = 0, V 2 7 ^ 0) to 
show that the QH state is stabilized for the smallest val¬ 
ues of V 2 as a result of the strain induced flat PLL-like 
bands. Second, we will look at the case (H 7 ^ 0, V 2 = 0) 
to show that a NN interaction will induce the ferroelec- 
tic charge ordered state. Third and last we will focus 
on selected cases of (Hi 7 ^ 0 , H 2 7 ^ 0 ) to show that the 
ferroelectric charge density wave is strongly suppressed 
as compared to the valley-ordered QH state, precisely 
due to the localization of the low-energy states flat band 
states on a single sublattice. 

Figure shows the HF results for various values of 
NNN interaction H 2 while keeping Hi = 0. Panels (b)-(c) 
show the QH order parameter Aqh( 0 defined in Eq. (28) 
as function of then tripled cell index 1. Black and red 
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FIG. 7. Panel (a) shows the maxima of QH order param¬ 
eter Aqh as function of NNN interaction H2 for various val¬ 
ues of the pseudo-magnetic field amplitude: A = 0.10 (red), 
A = 0.15 (blue), A — 0.20 (black). Panels (b)-(c) show the 
QH order parameter as function of tripled unit cell Aqh(1) 
(Z labeling tripled unit cell) for the same values of A [ (b) 
A = 0.10; (c) A = 0.15; (d) A = 0.20]. Curves are shown 
for H 2 = 0.1,0.2,0.3, 0.4, 0.5, 0.6, 0.7,0.8,0.9,1.0, in descend¬ 
ing order, i.e. bottom most curves H 2 = 1.0 and top most 
curves V 2 = 0.1. Black and red curves correspond to A and 
B sublattice, respectively. Hi = 0 in all cases. 


curves correspond to A and B sublattice, respectively, 
and the strength of the tunneling amplitude variation A 
is A = 0.10 (b), A = 0.15 (c), and A = 0.20 (c). It is 
apparent from these figures that the QH order param¬ 
eter follows the profile of the effective pseudo-magnetic 
field. In the spatial region where the flat band states 
are localized on the A sublattice, the QH order param¬ 
eter develops predominantly on the A sublattice, and 
vice versa for the B sublattice region. Ordinarily, in the 
honeycomb lattice QH state, the spontaneously induced 
magnetic fluxes are opposite on the A and B sublattice, 
averaging to zero over an elementary unit cell [7]. In 
the present case the localization of flat band states on a 
single sublattice leads to finite fluxes in regions of pos¬ 
itive and negative pseudo-magnetic field, which average 
to zero only over the larger strain superlattice unit cell. 

In addition to the locking of the QH order parameter 
Aqh( 0 fo the sublattice structure of the flat band zero 
modes, we observe that the strain-induced reorganization 
of the low-energy electronic structure into PLLs, funda¬ 
mentally changes the impact of interactions. Figure]^ a) 
shows the dependence of the QH order AQH(^max); where 
^max is the cell index where the QH order is strongest, on 
the strength of the interaction H 2 for various values of 
the pseudo-magnetic field strength (i.e., A). Whereas 
for the unstrained honeycomb lattice a finite interaction 
H 2 ~ 1.3 is needed to stabilize the QH state (within HF 
theory) [43H45] due to the vanishing density of states at 
half filling, here we find that the QH state is induced 
already for small interactions, in particular in regions 
where the effective field is strongest. For perfectly flat 
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FIG. 8. (a) Plots of the charge density wave order pa¬ 

rameter pa{1) — Pb{1) (red) and total charge redistribution 
Pa{1) + Pb{1) — 1 (black) for values of the NN interaction 
Vi = 0.2,0.4,0.5, 0.6, 0.7, 0.8, in both cases plotted in ascend¬ 
ing order (top curve Vi = 0.8, bottom Vi = 0.2); A = 0.10. 
(b) Same is (a) but with A = 0.20. 


bands such as PLLs one expects the interaction-induced 
order to scale linearly with interaction for weak coupling 
[iniHT]. This is reflected in Fig.j^a) which shows that for 
increasing pseudo-magnetic field, setting the PLL band 
flatness, the QH is more robust and its dependence on 
the interaction is approximately linear, with deviations 
at very small values of V 2 ■ 

Next, we turn to the case of finite NN interaction Vi 
while keeping V 2 = 0- In the absence of strain the unfrus¬ 
trated NN interaction will favor a CDW characterized 
by translational symmetry preserving sublattice charge 
imbalance [43ll45] . On the contrary, in the presence of 
strain, adopting the PLL picture for the low-energy elec¬ 
trons and focusing on the zeroth PLL, the effect of the 
NN interaction is expected to be suppressed, as the ze¬ 
roth PLL states live exclusively on one sublattice. Nev¬ 
ertheless, since the NN has the potential to cause charge 
asymmetry between the sublattices (i.e., the regions of 
positive and negative pseudo-magnetic field), and higher 
PLLs may be relevant to the energetics, we anticipate 
the system to develop a charge density wave of ferroelec¬ 
tric type (i.e., Ap+), with excess charge in regions where 
the pseudo-magnetic field is positive, and defect charge 
where it is negatice, or vice versa. In addition, in the 
previous section we observed that under the assumption 
of a uniform CDW the strain-induced PLL spectrum is 
gapped out. 

In Fig.j^we show results for both Ap+(Z) and Ap_{l), 
defined in Eqs. ( [^ and ( [2^ , for different values of the in¬ 
teraction Vi and strain A. As expected, we find the CDW 
order parameter Ap_(l) (shown in red) to become finite 
in regions where the pseudo-magnetic field is strongest, 
but to have the same sign in both positive and negative 
regions. The sign of concomittant ferroelectric polariza¬ 
tion depends on the sign of the pseudo-magnetic field. In 
the region where the field is positive the flat band states 
are localized on the A sublattice and pushing them down 
in energy, signaled by positive Ap_(?), leads to excess 
charge in that region at the expense of charge in the re¬ 
gion of negative field. At the same time we observe that 
the FP is more pronounced for stronger strains, and that 



FIG. 9. Plot of both the QH order parameter Aqh(0 (only 
A-sublattice shown; red dot curves) and total charge redis¬ 
tribution pa{1) + Pb{1) — 1 (black dot curves) as function of 
tripled unit cell index n for A — 0.20 and V 2 = 0.20. Different 
curves correpond to different values of Vi, explicitly labeled 
for clarity. 

it is very weak for small interaction. The latter may be 
attributed to the fact that NN interactions have no effect 
in the zeroth PLL. 

We proceed to consider the case of both finite Vi and 
V 2 - We have seen that both of these interaction individ¬ 
ually favor different gapped ground states. At the same 
time we argued that as a consequence of the different na¬ 
ture of these interactions, i.e. Vi being inter-sublattice 
and V 2 being intra-sublattice, they have different impact 
on the low-energy flat band electrons. Due to the spatial 
separation of states localized on different sublattices, it is 
expected that the effect of Vi is suppressed. One there¬ 
fore expects the QH state to survive even for V 2 < Vi, 
which corresponds to the physically relevant regime. 

In Fig. we present results for various Vi with finite 
V 2 = 0.2 and for A = 0.20. We plot both the ferroelec¬ 
tric order parameter Ap_|_(/) and the QH order param¬ 
eter Aqh( 0 for the A sublattice. The key observation 
is that even for small V 2 = 0.2 the QH survives up to 
NN interactions Vi ~ 1.0. The therefore conclude that 
the effect of interactions on periodic strain induced flat 
bands follows from their PLL character. In particular, 
the sublattice structure of the low-energy flat bands is 
the decisive factor in determining which order is sponta¬ 
neously generated by interactions. 


B. Remarks on spinful electrons 

We close this section with a number of remarks on the 
electron spin. The numerical calculations have been per¬ 
formed for spinless electrons in graphene. Taking spin 
into account gives rise to a richer structure of polarized 
states, specifically the ferromagnetically (FM) and anti- 
ferromagnetically (AFM) polarized states should then be 
considered, in addition to the Quantum Spin Hall polar¬ 
ized state. Moreover, the argument for the suppression of 
the NN interactions does not apply to the onsite Hubbard 
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interaction, Hu = U which must be included 

in the interacting Hamiltonian. 

The results for spinless electrons in graphene do, how¬ 
ever, directly apply to TCI surface, which do not have an 
additional degenerate spin degree of freedom. Instead, 
as a consequence of spin-orbit coupling, spin is already 
part of the low-energy Dirac structure. In particular, 
our results imply that on the surface of a TCI and in the 
presence of periodic strain, interactions will lead to the 
formation of the QH state. 

As a first step towards understanding the polarization 
of spinful strain superlattice-induced PLLs in graphene, 
we have performed numerical Hartree-Fock calculations 
with an interacting Hamiltonian given by Hu- We 
find that the mean-field ground is a superlattice anti- 
ferromagnet. The superlattice anti-ferromagnet exhibits 
anti-ferromagnetic order, as expected on a bipartite hon¬ 
eycomb lattice, however, since the flat band states are 
localized on one sublattice only, in each of the two re¬ 
gions of the strain superlattice an effective magnetiza¬ 
tion develops. Hence, as a consequence of the particular 
structure of the zeroth PLL states, the anti-ferromagnetic 
order is transferred to the superlattice. 

A similar result was reported in Ref. [551 which found 
anti-ferromagnetic order induced by non-periodic strain, 
where the bulk and the sample boundary have an effective 
but opposite magnetization. 


VI. DISCUSSION AND CONCLUSION 

We have shown that in the presence of a strain su¬ 
perlattice, a periodic modulation of elastic lattice defor¬ 
mations, a system of low-energy Dirac electrons exhibits 
a fourfold degenerate zero energy flat band, reminiscent 
of a zeroth PLL. The PLL structure originates from the 
pseudo-magnetic field, generated by nonuniform strain. 
The strain superlattice unit cell consists of two spatially 
distinct regions, one in which electrons see a positive 
pseudo-magnetic field and one in which they see a nega¬ 
tive field. The single-particle states of the degenerate flat 
band have a special and important localization property: 
in each of the two regions their wavefunction has support 
only on one of the Dirac pseudospin species. 

Periodic pseudo-magnetic fields can occur both in 
graphene and on the surface of a TCI, which hosts 
pairs of Dirac fermions at opposite momenta related 
by time-reversal symmetry. The important fact that 
Dirac fermions are unpinned to time-reversal-invariant 
momenta in the BZ allow for pseudo-gauge field under 
time-reversal invariant perturbations such as strain. 

Interactions between electrons in continuum PLLs are 
expected to lead to the formation of polarized states split¬ 
ting the degeneracies of PLLs. We have investigated PLL 
polarization for the case of lattice PLLs corresponding 
to periodic pseudo-magnetic fields. Two polarized states 
were shown to fully lift the zero energy flat band de¬ 
generacy while at the same time pushing all occupied 


(unoccupied) lattice PLLs down (up). The first state is 
the sublattice polarized charge-ordered state, which can 
be pictured as a spatially polarized state with all PLLs 
of the positive (or negative) field region occupied. The 
second is the anti-ferro-valley ordered state, or sponta¬ 
neous quantum Hall state, for which all PLL states effec¬ 
tively seeing a positive (or negative) field are occupied, 
implying time-reversal symmetry breaking. We found 
that other polarized states, even though they have simi¬ 
lar characteristics in the continuum, do not fully lift the 
lattice flat band PLL degeneracies. 

Using self-consistent Hartree-Fock calculations we have 
studied an interacting honeycom lattice model with pe¬ 
riodic strain-induced lattice PLLs. Our results demon¬ 
strate that the strain-induced reconstruction of low- 
energy elecronic structure, in particular the presence of a 
zero energy flat band, determines the impact of interac¬ 
tions. Three key results highlight this conclusion. First, 
the mean-field order parameters clearly reflect the pe¬ 
riodicity of the pseudo-magnetic field, showing that the 
amplitude of the order parameter is tied to the strength 
of pseudo-magnetic field. 

Second, as a consequence of the characteristic wave- 
function support of the flat band single-particle states, 
the effect of interactions that favor time-reversal invari¬ 
ant pseudo-spin order is suppressed. This effectively en¬ 
hances interactions that favor time-reversal symmetry 
breaking valley order with associated spontaneous quan¬ 
tum Hall effect. We have established this result in the 
context of the graphene lattice model, where the pseudo¬ 
spin corresponds to the sublattice degree of freedom. 
The NN interactions are inter-sublattice interactions, and 
therefore supressed, as the flat band single-particle states 
have support on one of the sublattices only, in each pos¬ 
itive or negative field region. The result, however, is 
general and applies equally to TCI surface states. The 
physical interpretation of pseudo-spin and valley are dif¬ 
ferent, as explained in Sec. [Hj yet the effect of interac¬ 
tions favoring time-reversal symmetry breaking remains 
strongly enhanced. 

Third, the valley-ordered spontaneous quantum Hall 
state already occurs for small interactions when the PLLs 
are well-developed. The PLL structure is a way to signif¬ 
icantly enhance density of states near the charge neutral 
point. We conclude that in the presence of periodic strain 
and interactions, a system of unpinned Dirac electrons 
has a generic instability towards a spontaneous quantum 
Hall phase. 

The emphasis of the numerical calculations we report, 
has been on spinless electrons in graphene. In TCI sur¬ 
face states, however, no additional spin degeneracy is 
present. Due to spin-orbit coupling the electron spin is 
an intrisinc part of the low-energy Dirac structure. In 
particular, this implies that there are no purely spin- 
polarized phases, such as the global anti-ferromagnet, 
competing with the spontaneous quantum Hall phase, 
favoring the latter as ground state. 

Whereas uniform strain-induced pseudo-magnetic 
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fields suffer from implementation limitations, particu¬ 
larly beyond the nano-scale, periodic strain can poten¬ 
tially be realized in macroscopic sample sizes. In fact, 
such periodic strain fields and induced pseudo-magnetic 
fields were demonstrated in TCI heterostructure [5^155] . 
making TCI surface states the prime candidate to exhibit 
spontaneous formation of nontrivial electronic states. 
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Appendix A: Landau levels for Dirac electrons 

For the purpose of being selfcontained we collect some 
standard results of Dirac fermions in a constant magnetic 
field in this Appendix. These may be directly applied to 
the case of time-reversal invariant pseudo-magnetic fields, 
bearing in mind the key characteristic of opposite sign of 
the pseudo-magnetic field in the two valleys. 

In the presence of a magnetic field we define the dy¬ 
namical momenta using the Peierls substitution 

Pa -t fla=Pa- eAa{r) = -iKda + |e|Ac(f), 

where Pa is the momentum operator and a = x,y (we 
restrict the description to two dimensions), fa are the 
position operators obeying [fa,Pi 3 \ = ih and Aa{f) is 
the electromagnetic gauge field. In a magnetic field the 
momentum operators IIq do not commute but instead 
obey the canonical commutation relation 

[na,n;3] = \e\[pa,Ai3{f)] + \e\[Aa{f),py] 

= -ih\e\Fap 

where Fay = daAy — dyAa is the field strength. Assum¬ 
ing a uniform field strength in the z direction the mag¬ 
netic field is given by B\ = exfivF^,^/2, implying that 
Ffj,i, = CfiiyxBx. In particular we have for a uniform field 
B = Bz in the z direction 

[FIq,,!!/?] = -ih\e\sgn{B)Beayz- (AI) 

In this expression we have explicitly separated the sign 
of the magnetic field from its strength B = \B\ so as to 
make dependencies on the sign of the field transparent. 
Defining the fundamental characteristic length scale in 
the system, the magnetic length, as Is = x/h/{\e\B), we 
can write [IIa,II^] = —ih'^sgn{B)€ayz/lB- This implies a 
canonical commutation relation between the dynamical 
momenta and inspires to define creation and annihilation 


operators in the usual way as 

+ fsgn(R)ny), 

a = - *sgn(R)ny), (A2) 

which obey [a, 6,^] = 1. Note that the definition of these 
operators depends on the sign of the R-field, which is a 
direct consequence of Eq. ( |Al[ ). We note in passing that 
all of the above did not require specifying a gauge for Aa ■ 
The Landau level spectrum of a Dirac Hamiltonian of 
the form 


B xQx F 


is then straightforwardly obtained by making the substi¬ 
tution hQfj, —>■ n^. Squaring the Hamiltonian yields 

= vUfil + H^) + 4[n„ Hj^jF^r^ 

= - I -f sgn{B)Tz) 

‘b 

We use that d^d = n for standard oscillator wave func¬ 
tions ifn, i-e. d^dpn = npn, dipn = x/npn -1 and 
d^Pn = x/nF Ipn+i- One therefore obtains the Landau 
level energies 


E±(n) = Fx/^, n = l,2,..., 
2 _ 


(A3) 


Each of the E± (n) is two-fold degenerate because of the 
valley degree of freedom, in addition to an Ny = Aj27rl^ 
degeneracy where A is the area of the system. 

We find the corresponding eigenstates by taking a 
closer look at the explicit expression for the Hamiltonian 
in each valley. Writing for the Hamiltonian in valley 
ly = F we have 


n. = iy 




The eigenstates belonging to the eigenvalues E±(n) of 
Eq. @ are easily obtained as 




1 / Iv^n—l,fc) \ 

\/2 \ ±v\pn,k) J 


(A4) 


In addition to the states |'I'„;,±) (n = 1,2,...) there 
are also zero mode states l^^oiy), one for each valley, which 
have zero energy (Eq = 0). Inspecting of Eq. (18) reveals 
that these states are given in each of the valleys as 





(A5) 
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We stress that this implies the zero modes are localized 
on opposite sublattices for the two valleys, as we had 
exchanged sublattices for the K- valley. 

We proceed to consider the effect of symmetry breaking 
terms on the Landau level spectrum. Specifically, we 
consider first the set of time-reversal symmetry invariant 
masses rh which enter the Hamiltonian as 

= TO • f = miiy^ + m2V^ 

For small masses we may use perturbation theory to 
study the splitting or shifting of Landau levels. It turns 
out however that the exact energies can be obtained in 
the presence of Hm- The energies are found by squaring 
the Hamiltonian % = Ho + which gives 

= Vp(fi^ + fly) + Vp\fix, ny]r2,ry -i- 
= - 1 -I- sgn(H)rz) -|- to^, 

where we use the anticommutation relations of the F- 
matrices. We directly find the energies 

E±{n) = ±^/2^'^n + m?, n = l,2,... (A6) 

Expanding the square root \j2^n^\ m? 12^'^n in 
small w?/2^^n yields the same result as second order 
perturbation theory. 

The Landau level spectrum in the presence of masses 
can alternatively obtained by using the relation [H*, Fj] = 
2ieijkl^k to construct a unitary matrix U which rotates 
the vector to so that one has U^m - YU = toFs, with 
TO = |to|. Such transformation block diagonalizes the 
Hamiltonian and we obtain the energies in a more direct 
way. In particular, we can employ the unitary rotation to 
find the energies of the zero modes for the case of massive 
Dirac fermions. In the rotated basis we can construct 
zero mode states in the same way as before, which will 
have energies 

Eqi, = -vm (A7) 

where v = ± represents the valley degree of freedom. 

We conclude by taking into account the time-reversal 
symmetry breaking but chiral symmetry preserving mass 
r] entering as Hjj = rjT^. Since such this term is a scalar 
under chiral rotations generated by H* we may consider 
H = T-Lo + 'Hfh + Hrf and use the chiral rotation U to 
block diagonalize the Hamiltonian. The total mass term 
then is mv^r^ which directly leads to the energies 

Ev±{n) = -k\/2^^n + {vm + rf)'^ (A8) 


Appendix B: Setup of Hartree-Fock calculations 


Our Hartree-Fock calculations in the presence of the 
strain superlattice and with interacting Hamiltonian (21) 
follow the scheme of Ref. HU In particular, we choose the 
unit cell of the unstrained lattice such that it contains 
six honeycomb lattice sites, as show in Fig. [^of the main 
text. This allows intra-sublattice mean-field structures 
to form, corresponding to modulation vectors K±, which 
connect the Dirac points of the honeycomb lattice. In 
terms of the elementary graphene lattice vectors 

di = a(l,y2)/2, 02 =a(l,-A/2)/2, 


the lattice vectors of the six-site unit cell are given by 
bi = 2di -I- 02, 62 = di — 02 . 


They are shown in Fig. together with the folded BZ. 
Figure also shows how, in the presence of the periodic 
strain superlattice, the superlattice unit cell is defined. 
The superlattice vectors are 61 and A 62 , in terms of the 
superlattive wave length A. 

The electronic Hamiltonian is expressed in terms of 
the fermion annihilation (and corresponding creation) 
operators Here a labels th sublattice {A/B), 

i = 1,2,3 labels the three sites of each sublattice species 
in the six-site unit cell, and I = 1,..., A lables the (six- 
site) unit cells in the superlattice unit cell, and ^ is a 
position index for the superlattice unit cell. The Fourier 
transforms is defined as 


The interacting Hamiltonian of Eq. (21) consists of two 
terms: the NN interaction and the NNN interaction. In 
momentum space the NN interaction Hy^ is given by 


kk' q 

+ (Bl) 

where repeated indices are summed. The matrix function 
X{q) connects NNs, and it is convenient to decompose it 
in the following way 

+ X°iq)Si,i, + X+{q)6pi,+y 


Clearly, X°j{q) connects sites within the same (six-site) 
unit cell, and it is explicitly given by 


for the Landau levels n = 1,2,.... The presense of both 
of these masses leads to a splitting of Landau levels, 
whereas the presence of either only shifts the energies. 
The presence of a time-reversal breaking mass breaks 
particle-hole symmetry in the n = 0 Landau level. Specif¬ 
ically, the zero mode energies become 

(A9) 


/ 1 0 1 \ 

X°{q) =1 10. (B2) 

1 1 / 

The functions X)^{tf) connect sites in different (six-site) 
unit cells and each only have a single nonzero entry. They 
are given by X+{^ = e-dH+b 2)-9 and X^^i^ = A 


Equ = -vm - 77 . 
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Similarly, the NNN interaction Hamiltonian is given 

by 


q)yi,,{q) 

kk'q 

+ q) (B3) 


The functions y‘^{q) connect NNNs on each sublattice 
a. They are directly obtained from Ref. HU taking into 
account the additional superlattice index I (in the same 
way as for X{q)). 

The quartic interactions of the interacting Hamiltoni¬ 
ans, schematically written as are decoupled in 

the standard mean-field way as (written schematically) 




the first line representing charge density order and the 
second bond density order. 

In terms of the actual superlattice electrons, the charge 


density order parameter is defined as 


1 


Pa^{l) = ^ k)). (B4) 


Bond order parameters are defined as straightforward 
generalizations of Ref. HD Of particular interest in our 
case is the QH order parameter, constructed from NNN 
bond order and defined by Eqs. (|2^ and (28). To obtain 


Eq. (27) we decompose y^ ji>{k — k'), wich arises due to 
the reordering of (B3), as 


yiA^-k') = Y. QAA)QZA-k'). ( b 5 ) 

/i=i 

The sum over /i follows from the three ways in which 
NNN sites may be connected. Note that fc) = 

Qttji'ik)- We can now NNN bond order parameters used 
to construct the QH order parameter of Eq. (28). The 
NNN bond order mean-fields are defined by 

A31' = ^ E k)A{i\ A- (B 6 ) 


The QH order paramater Aqh( 0 is defined so that each 
unit cell labeled by I is associated with 2x9 NNN bonds. 
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